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Abstract. Wc present the first systematic derivation of tlic one-loop correction to 
the large scale matter power spectrum in a mixed cold+hot dark matter cosmology 
with subdominant massive neutrino hot dark matter. Starting with the equations of 
motion for the density and velocity fields, we derive perturbative solutions to these 
quantities and construct recursion relations for the interaction kernels, noting and 
justifying all approximations along the way. We find interaction kernels similar to 
those for a cold dark matter-only universe, but with additional dependences on the 
neutrino energy density fraction fi, and the linear growth functions of the incoming 
wavevectors. Compared with the /ly = case, the one-loop corrected matter power 
spectrum for a mixed dark matter cosmology exhibits a decrease in small scale power 
exceeding the canonical ^ 8/^ suppression predicted by linear theory, a feature also 
seen in multi-component A^-body simulations. 
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1. Introduction 

Standard big bang tlieory predicts a background of relic neutrinos, permeating the 
universe at an average of 112 neutrinos per cubic cm per neutrino flavour. This 
enormous abundance means that even a sub-eV to eV neutrino mass m^, will render these 
otherwise elusive particles a significant dark matter component, Vtyh"^ = J^f^u/i^'i eV). 
Experimentally, a lower limit of J2^u ^ 0.05 eV on the sum of the neutrino masses 
has been established by neutrino oscillation experiments (e.g., [1]). On the other 
hand, tritium /3-decay end-point spectrum measurements point to an upper limit of 
YlfT^u ;$ 6 eV (e.g., [1]). The corresponding energy densities, 0.001 < f^j^ < 0.12, make 
for an interesting prediction to be tested against cosmological observations. 

Importantly, neutrino dark matter is of the hot variety because of their inherent 
thermal velocity, which prevents them from clustering gravitationally on scales smaller 
than their free-streaming length. This free-streaming effect then feedbacks into the 
evolution of the dominant cold dark matter (CDM) component via the gravitational 
source term, leading to a suppressed rate in the formation of structures on small length 
scales that is in principle manifest in the power spectrum of the large scale structure 
distribution [2-6]. At present, observations of the cosmic microwave background 
anisotropies, galaxy clustering, and type la supernovae together engender a conservative 
upper limit on the contribution of neutrino dark matter, or equivalently, on the neutrino 
masses, of < 1 eV within the ACDM framework. See [7, 8] for recent reviews. 

Many future cosmological probes will continue to improve on this limit, or perhaps 
even detect neutrino dark matter (e.g., [7,9-17]). To this end, it is important to 
note that many of these probes, particularly high redshift galaxy surveys and weak 
gravitational lensing, will derive most of their constraining power at wavenumbers 
k 0.1 1 h Mpc~^, where the evolution of density perturbations has become weakly 
nonlinear. Coincidentally, these are also the scales at which neutrino free-streaming is 
expected to produce the largest effect on the large scale matter power spectrum. Thus, 
it is crucial that we understand the nonlinear evolution of density perturbations at these 
scales, in order to maximise our gain in detecting/constraining neutrino dark matter. 

A number of recent works have attempted to model the effects of massive 
neutrinos on the matter power spectrum at nonlinear scales using various different 
techniques. These include multi-component A^-body simulations [18] and semi-analytic 
halo models [11,19,20]. However, with the exception of A^-body simulations, these 
methods all contain elements that require prior calibration against simulation results, 
and not surprisingly, all calibrations to date have been performed in a ACDM setting. 
This renders these nonlinear models at present not completely satisfactory for use with 
cosmologies containing massive neutrinos. Recalibration against simulations in the 
appropriate cosmology will, however, enhance/restore our confidence in these methods. 

Recently, Saito et al. [21] proposed an alternative method based on higher order 
cosmological perturbation theory. Since its inception in the 1980s [22-25], higher 
order cosmological perturbation theory has found numerous applications ranging from 
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computation of the weakly nonlinear power spectrum and gravity-induced bispectrum, 
to the exploration of nonlinear galaxy bias. See reference [26] for a review. In the 
perturbative approach, one envisages nonlinear evolution as the outcome of interactions 
of a collection of linear waves (perturbations). The equations of motion of the system 
define the interaction kernels. 

Saito et al.^s recipe is simple: assume the neutrino density perturbations remain 
linear at all times, and apply nonlinear modelling only to the CDM+baryon component. 
This basic scheme has also been adopted in some earlier nonlinear models including 
massive neutrinos [11]. For J2^u ^ 0.6 eV, multi-component A^-body simulations have 
confirmed that a linear approximation for the neutrino density contrast is sound [18]. 

For the CDM+baryon component, Saito et al. calculated the one- loop correction 
to the CDM+baryon auto-correlation power spectrum, using the correctly computed 
linear waves (perturbations), but interaction kernels that have been developed for a 
CDM-only universe. This amounts to ignoring additional mode-coupling effects between 
the linear growth functions at different wavenumbers. Although, as we shall show, this 
approximation does lead to a considerable simplification in the final form of the nonlinear 
power spectrum, it is nonetheless reminiscent of the mismatch between the cosmology 
to which we apply and that on which we calibrate some of the semi-analytic methods 
discussed above, and therefore calls for closer scrutiny. 

In this paper we present a systematic derivation of the one-loop correction 
to the matter power spectrum in the presence of massive neutrinos from first 
principles. As in [11,21], we assume linearity for the neutrino component. For 
CDM+baryons, however, we begin with the relevant equations of motion, and solve 
them (approximately) perturbatively in an Einstein-de Sitter VLm = 1 background 
at high redshifts [z > 0.5). From these solutions we derive interaction kernels that 
capture also the physics of mode-coupling between linear growth functions at different 
wavenumbers. With these kernels we compute the correct nonlinear power spectrum. 

In section 2 we give the relevant equations of motion. We discuss briefly the linear 
order solutions in section 3, before formulating our higher order perturbative approach 
in section 4. In section 5 we outline our scheme to obtain approximate solutions to the 
equations of motion for higher order perturbations, while in section 6 we derive recursion 
relations for the interaction kernels and thus generalise our approximate solutions to 
arbitrary orders in perturbative expansion. Section 7 deals with the calculation of the 
one-loop correction to the matter power spectrum, which we evaluate numerically for 
realistic cosmological models and discuss in detail in section 8. We conclude in section 9. 

2. Equations of motion 

We begin with the standard set of equations of motion for the density perturbations 
S{x,t) and the peculiar velocity u{x,t) of the CDM+baryon component (e.g., [26]), 




+ V ■ {[l + 6{x,r)]u{x,T)} = 0, 



dr 
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^^^^^ + niT)u{x, t) + [uix, r) ■ S7]u{x, r) + V<f (a;, r) = 0, (2.1) 

where r is the conformal time, x the comoving coordinates, Ti. = dlna/dr = Ha the 
conformal expansion rate, and the stress tensor has been set exphcitly to zero. The 
Poisson equation 

V'^x,r) = ^n\r)Q.mirm - fu)Six,r) + fj^{x,r)] (2.2) 

relates the Newtonian gravitational potential $(a;,r) to the density perturbations. 
Perturbations from both CDM+baryons 6{x, r) and neutrinos S^lx, r) contribute to this 
gravitational source term, although in the latter case the contribution is suppressed by 
the neutrino density fraction fi, = r2,^(r)/f2m(r) = f2i,^o/^m,o- This fraction is constant 
in time once the neutrinos have become nonrelativistic. 

It is common to rewrite the equations of motion in Fourier space. Defining the 
Fourier transform 

~ f d^x 

(p{k, r) = J exp(-ifc ■ x)(l){x, r) (2.3) 

for some field (f){x, r), and the divergence of the velocity field 

e{x,T)=V ■u{x,t), (2.4) 
equations (2.1) and (2.2) can be equivalently expressed as 

d6{k,T) 



dr 



e{k, t) = - J d^qid^q26Dik - qi2)a(qi, q2)0{qi, T)6{q2, r). 



where 



+ n{T)9{k, r) + ^n\T)n„,{T) [(l - U)~6{k, r) + fj-^{k, r)] = 
- J d^q^d^q^Soi^k - qu)l3iqi, 92)^(^1, r)e{q2, r), (2.5) 

«(qi, 92) = PiQi, Q2) = ^00 • (2.6) 

Here, we have defined q = \q\, and qi...j = qi + ■■■ + qj. The function 6d denotes 
the Dirac delta function. As usual, we assume the vorticity of the velocity field 
w{x, r) = V X u{x, r) to be negligible at all times. 

For the neutrino component, a full treatment requires that we track the evolution of 
the neutrino phase space density /^^ (a;, p, r) by solving the Vlasov equation (e.g., [26]), 

^ + ^ ■ Vr - am,V<l> ■ ^ = 0, (2.7) 
or my a op 

where my is the neutrino mass, and we have assumed the neutrinos to be nonrelativistic. 
Integrating over the neutrino momentum p, 

j d'p r{x,p, r) ^ p^(r)[l + 5^{x, r)], (2.8) 

yields the neutrino density contrast 5'^(x,r). 
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3. Linear theory 

Linearising the equations of motion (2.5), i.e., dropping all terms on the r.h.s. of the equal 
sign, the growing mode solutions for the CDM+baryon density contrast and peculiar 
velocity at some wavevector k can be written as 

5i(fc,r) = Di(A;,r)5(fc,ro), 

^i(fc, r) = -n{T)f{k, T)Di{k, T)6{k, To), (3.1) 

where tq denotes some initial time well in the matter-domination epoch, Di{k, r) is the 
linear growth function, its logarithmic derivative 

^_ d\nD,{k,T) 1 d\nD^{k,T) 

= d\na = H 97^' ^^-2^ 

and we have dropped the tildes on 5i(fc,r) and ^i(fc, r) for convenience. A similar 
expression, 

5'ak,T) = D'ak,T)5''ik,To), (3.3) 

describes the growth of the neutrino density perturbations in the linear regime. 

The inclusion of massive neutrinos in the matter content of the universe introduces 
a new length scale to the problem, the neutrino free-streaming scale /cps (e-g-? [27]), 



where 



3C(3) _ 81 /eV\ 



\ , / N — / N 1 I km s . (3.5) 
\\ 21n(2) a(r) \mj ^ ' 

At wavenumbers k -C /cfs, neutrinos cluster gravitationally and behave essentially like 
CDM. In the other limit k ^ /cpg, the neutrinos' inherent thermal velocity prevents 
efficient infall into gravitational potential wells, thereby suppressing the neutrino density 
perturbations relatively to their CDM+baryon counterparts. 

This qualitative picture is generally true for any cosmology. However, if we restrict 
our considerations to an Einstein-de Sitter universe, then using the relations a oc 
and 7i(r) = 2/r it is easy to show that 

Di{k,T)(xa, /(/c,r) = l, k <^ kps, 

Di{k,T) oca^-", /(A;, r) = 1 - /i. A; > A;fs, (3.6) 

with 



The k ^ kps solution is obtained by setting (5^(fc, r) explicitly to zero in the equations 
of motion. For the neutrino component, the expressions 

D^,{k, t) D^{k, r) k « A^Fs, 

[K, To) 



6-{k,To) 



kj.sir){l-U) 

k"^ — kpa,{T)fi, 



k > A;fs (3.8) 
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have been shown to be asymptotic solutions to equation (2.7) [27]. For intermediate 
k values, the linear growth functions for both the CDM+baryon and the neutrino 
components must be calculated numerically with a Boltzmann code such as CAMB [28], 
which also gives the neutrino sector a full general relativistic treatment [29]. 



3.1. Fitting formulae 

We shall need estimates of the linear growth function and particularly its logarithmic 
derivative later in the analysis. For this we resort to fitting formulae. 
For the CDM+baryon component, we have [30] 

D,{k,r) ~ \l+r%k,T)f''D\->'{r), (3.9) 

and hence 

/(A;,r)~l-- , (3.10) 

where c ^ 0.7, and 



r{k, r) 



1+Z/Fs(fc)' 



y^s{k) = 17.2/,(1 + 0.488/-^/^)(piV,//,)2, 

P= {^^j^^Ql,{^^,oh'r\ (3.11) 

Here, Di (r) is the growth function in the absence of massive neutrinos, normalised such 
that -Di(t) = a/ceq in an Einstein-de Sitter universe. The quantity N^, is the number 
of massive neutrinos (assuming equal masses) which we generally take to be three, and 
62.7 is defined through Tcmb = 2.7 62.7 K. The claimed accuracy of the formula (3.9) 
is 1 to 2 % [30]. 

For the neutrino component, we find 



6-{k,T^) 



kUr){l-U 



[k + k^s{r)Y-kUT)U 



(3.12) 



to be a good interpolation between the k <^ A;fs and k ^ /cfs limits, accurate to better 
than 5 %. 



3.2. Linear power spectrum 

The linear matter power spectrum is defined as 

P\k,T)5D{k + k') = {6j{k,r)6^{k',T)), (3.13) 

where 6j{k,T) = fcbSi{k,T) + fu6i{k,T), with fcb = {^c + ^b)/^m, counts both the 
CDM+baryon and the neutrino density contrast. In terms of the linear growth functions, 

P\k, r) = f,,P^,{k, r) + 2/,,/,P,i(A;, r) + /JP,^(fc, r), (3.14) 
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with 

Pi(fc,r) = [D,{k,T)T,,{k,T,)fP\k), 

P^{k,T) = [D'i{k,T)%{k,T,)fP\k,), 

P^k, r) = D,ik, T)D'i{k, T)T,,{k, To)mk, r^)P\k), (3.15) 

where Ti{k,To) are the hnear transfer functions mapping the initial fluctuations (from, 
e.g., inflation) through the epochs of horizon crossing and matter-radiation equality to 
time To, i.e., 6{k, tq) = T{k, ro)5^(fc). We have also assumed in equation (3.15) adiabatic 
initial conditions so that P^bik) = Py{k) = P^f,^{k) = P\k). Note that the power 
spectrum defined in this manner contributes 47rk^P{k) per logarithmic wavenumber 
to the variance, in contrast to, e.g., the default output of CAMB, which contributes 
k^P{k) I (27r^) between In k and In /c + c? In /c. 



4. Beyond linear theory 



We wish to find a higher order perturbative description for the CDM+baryon density 
contrast and peculiar velocity. To do so it is convenient to define a new time variable 



s = In a(r) = 7-^(r)(ir, 
and the vectors 



'^^{k,s) 







n\s)^(^^ s) 

The equations of motion (2.5) can then be written in a more compact form, 

ds'^aik, S) + /C,^fc(fc, S) + Nab'^lik, S) = 



(4.1) 
(4.2) 



(4.3) 



where a,b,c = 1,2, and repeated indices imply summation. 

In the homogeneous part of equation (4.3), the matrix K is given by 

-1 



K 



(4.4) 



where we have assumed explicitly an Einstein-de Sitter universe. The inhomogeneous 
part is specified by the matrix, 





N 



(4.5) 



and the tensor 7abc(gi, ^2) is zero except for 7121 = a{qi, ^2) and 7222 = I3{q\, 92) 
We seek a perturbative solution of the form 

00 

^>{k,s) = Y.^^''\k,s). 

n=l 

For the neutrinos, we assume only the linear solution is nonzero, i.e., 

5l{k,s) 



*'^(fc,s) 







(4.6) 



(4.7) 
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Then the equation of motion for the nth order solution, where n > 2, is 

s) + K,,^l^t\k. s) = Bjr\k, s), (4.8) 

with 

B^J'\k,s)= J d^qid^q25Dik - qi2habciqi,q2) 

n-l 

xY.i^t"'\qi,s)4"'\q2,s). (4.9) 

m=l 

Observe how the term proportional to '^'^{k,s) has disappeared for n > 1 because of 
the assumption that neutrino density perturbations remain linear at all times. 



5. Approximate solutions 

Equation (4.8) can be solved by first defining a transformation matrix U via 

K = U-^KU = Diag(Ki, K2), (5.1) 

so that in the diagonal basis, the equation of motion becomes 

d,^lr\k, s) = -K.#)(fc, s) + Bj^\k, s), (5.2) 

where ^/''•"^ = U^^il)^"'\ and 5^"^ = U~^B^^\ The formal solution to equation (5.2) is 
simple, 

i^^J'^k, s) = e'^^o ""^'Vi"^^, So) + r e-T' B^J'\k, s')ds' 

J so 

= r e-I>'^'''"Bi"\k,s')ds', (5.3) 

J so 

where we have used 'ip^\k,So) = in the second line, i.e., the perturbations at the 
initial time Sq are completely described by linear theory. Within our formulation of the 
problem the solution (5.3) is exact. To gain further ground, however, we must make 
some approximations, as we show order by order below. 

(i) The n = 1 solution is by definition 

^W(fc,s) = J d'q6n{k-q)Ql^\q,s)6,{q,s), (5.4) 

where Q^i\q, s) = 1, and Q'i\q, s) = f{q, s). 

(ii) For n = 2, we have 

(fc,s)= J d^qid^q2SDik-qi2)Di{qi,s)Di{q2,s) 

X U~^^-fbcd{quq2)Qi^\qi,s)Qi\q2,s)6{qi,So)S{q2,So). (5.5) 
The solution in the diagonal basis is then 

= f'ds' I d^q,d^q2 e-/>'"^^"e[''(^^'^')+''('?^'^')l^' 

J so J 

xSoik - qi2)U-f,^-fbcd{qi, 92) 

xQW(gi, s')Qi\q2, s')5(qi, ^0)^(^2, Sq), (5.6) 
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where we have defined s) = \nDi{q, s)/ Ina(s). 

The integration over s' can be performed using integration by parts, i.e., 

ds'Q^^\q^, s')Qi\q2, s')e-Is' -ads" ^Mg„s')+H,„s')]s' ^ 

Q^cKqi,s')Qa\q2,s') ^_ £^ Kads" ^[h{qus')+h{q2,s')]s 

' ds'Q<^^\q,, s')Q^J\q2, s')e-Is' -^ds" ^[h(qus')+Hq„s')W 



so 



so 



d 



K.a + fiqi,s') + /(g2,s') 



(5.7) 



Q^^\qus')Q^\q2,s')ds' 

and we have used the relation h{q, s) + {dh/ds)s = f{q, s). Compared with the original 
integrand on the l.h.s., the integrand on the r.h.s. of equation (5.7) is suppressed by the 
time derivative of f{q,s), whose size we estimate using the fitting formulae (3.9) and 
(3.10) to be 

df{q,s) ^ c^ir' 

ds "(l + r-)2' ^''■''^ 

and Of /ds > for all s, meaning that the linear growth rate cannot decrease in 
an Einstein-de Sitter cosmology. For a given q value, df /ds peaks at r = 1, or 
a = [l+yFs(g)] CLeq, wheu the neutrinos transit from a non-clustering to a clustering dark 
matter, thereby enhancing the (linear) gravitational source term. Here, df /ds\r=i — 
cfi/4: ^ 0.1 fi^. Since the integrand is positive definite, we can conclude based on this 
estimate that the absolute fractional error incurred by dropping the second term of the 
integral (5.7) is bounded from above by ~ 0.1/j. (worst case: ~ —1, c = 1, d = 2, 
Qi = ^2); the actual error is likely to be much less. Hence we shall proceed by keeping 
only the first term of equation (5.7). 

Of what remains of the time integral (5.7), we retain as usual only the non-decaying 
part, leading to 

Di{qi,s)Di{q2,s)f{qi,s) 



d-^qid-^q2- 



iia + fiqi,s) + /(g2,s) 
xSoik - qi2)U~h^-fbcd{qi, qi) 

Q'-cHqi, s)<5d ^(92, s)S{qi, So)6{q2, Sq). 
Rotating back to the original basis gives us 



where 

Q^^\qi,q2;s)^ 
is the interaction kernel. 



J d^qid^q25D{k - gi2)Qi^'(9i, ^2; s)5i(qi, s)5i{q2, s) 

UabU^c^-fcdeiqi, q2)Qi\qi, s)Q'^\q2, s) 



Kb + f{qi,s) + f{q2,s) 



(5.9) 
(5.10) 

(5.11) 
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(iii) For n = 3, the same procedure and approximation scheme yield 

xQa\Qi^ 92, qs] s)6i{qi, s)6i{q2, s)6i{q3, s), 



with 



Qa\Qi, Q2, q-s] s) 



lib + /(gi, s) + /(g2, s) + /(ga, s) 



(5.12) 



X [7cde(9i2, q3)Qd\Qu 92; s)Q''^\q3, s) 

+ lcde{qi,q23;s)Q'a\qus)Qf\q2,q3;s)] . (5.13) 

The maximum absolute fractional error incurred by neglecting the time derivatives of 
/(g, s) is again estimated at a negligible ~ 0.1/;, for the n = 3 solution. 



6. General nth order solution and recursion relations 

Generalising to the nth order, we obtain the approximate solution 
ijj^\k, t) = J d'qi ■ ■ ■ d^qnSnik - q,...^) 

xQi^H^i, ■ ■ ■ , 9n; r)5i(qi, r) ■ ■ ■ 5i(q„, r), (6.1) 
where the interaction kernel Q*-"^ can be constructed from the recursion relation 



Qt\(li^ ■ ■ ■ , 9n; r) = a^fe^(gi, ■ ■ ■ , gn; r 



n-l 



'Jbcd(^qi---mi qm+l-- n)Q c (^l, ' ' ' 5 9m) t) 



m=l 



9m+l, 



9n; r j 



with 



<6 (gi, 



-gnj 



Kc + t^(")(gi,---,gn;r) 



2cu(") + 1 2 
3(1-/,) 2^W 



and 



Ar(")(gi, ■ ■ ■ , g„; r) = {2uj^^^ + 3)(c^(") - 1) + 3/,, 

" 91nA(g^,r) 



^^"^(gi,---, 



gn;rj 



E 



(91nafr) 



(6.2) 



(6.3) 



(6.4) 



m=l m=l 

Equations (6.1) to (6.4) should be compared with, e.g., equations (41) to (44), (83) and 
(84) of reference [26] for a CDM-only universe. 

In the limit /, — > 0, we have a;'^"'^(gi, ■ ■ ■ , gn; t) n, so that the a*^") matrix (6.3) 
depends only on n such as in equation (84) of [26]. This renders the interaction 
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kernel (6.2) into the standard expressions Fn{qi, ■ ■ ■ , q„,) and G„(qi, ■ ■ ■ , qn) for a CDM- 
only universe given in equations (43) and (44) of [26], i.e., 

Qt\Qi, Qn] t) Fn{qi, ■ ■ ■ , qn), 

#(9i, ■ ■ ■ , r) ^ Gn{qi, ■ ■ ■ , qn)- (6.5) 

We caution at this point that the terms proportional to in equations (6.3) and (6.4) 
count only the explicit dependences on f^; the function u^^^ also depends implicitly on 
through the linear growth function. 



6.1. Symmetrised kernels 



In practice we use the symmetrised kernels Q^"'\qi, ■ ■ ■ , q„,; r), constructed by summing 
Q^"'\qi, ■ ■ ■ ! 9n; t) over all permutations of the momenta gi, ■ ■ ■ , 9n and then dividing 
by n\. For future reference, we give here exphcit expressions for Q^^-*(qfi, 92; t) and 
Q^^'^{qi, qi q^] t), cast in a form as close to the standard CDM-only ones as possible: 



Qi {qi, q2, r) = -Ai + 

7 2 qiq2 



Q2 [Qi, 92; r) = -Ci + 

7 2 qiq2 



12 



(I2 



C2 \- C*3 — 

12 qi 



^2 {q, ■ q^)' 
^4 (gi ■ g2)' 

7 gig2 



(6.6) 



with 



^2 

C4 



^^i?(gi,g2)[/(gi) + /(g2)], 

/ (92) Wii (gi , ^2) + f^i? (gi , 92 ) / (gi ) ] , 
/(gi)[f^i?(gi> <i2) + c^i?(gi, g2)/(g2)], 

242^(^1' g2)/(gi)/(g2), 

^4?(gi,g2)[/(gi) + /(g2)], 

/(g2)[4?(gi:g2) + 4?(gi>g2)/(gi)], 
/(gi)[4?(gi,g2) + 4?(gi,g2)/(g2)], 

-4J(gi,g2)/(gi)/(g2), 



(6.7) 



where we have dropped the r label in cr^'^\qi,q2]'T) and f{q,T) for convenience. The 
factors v4i^2,3,4 and Ci, 2,3,4 have been defined so that they tend to unity as 0. 
The symmetrised form of Q*-^^ can be written as 



Q?{qi,q2,q3;r) 



+ 



7 



123 



.1 Ix 

3 \l8 

gi23 • qi2 



9123 ■ qs 



123 



Q?\qi,q2]r) 



4 



-z 



123 



9^231912 • 93) 
^12^3 



Q2\qi,q2]r) 



+ cyclic permutations \ 



(6.8) 
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where 

18 

Xi23 = X{qi,q2,q-i) = y c^n (^i, ^2, g3)/(g3), 

18 

Yi23 = Y{qi, gs, Qs) = yf^n (^i, ^2, gs), 

Zu3 = Z{qi, g2, ga) = Qf^uiQi, ^2, q3)f{q3)- (6.9) 
Again, X123, ^123, ^123 ^ 1 as /j, ^ 0. To obtain Q'i\qi, ^2, ^3; 1"), simply replace a[^^ 

(S) (3) (3) 

with (721 , and cr{2 with 0-22 the expressions for X123, ^123, and ^123. 



7. Power spectra 

We are interested in the total matter power spectrum P{k,T). As in the case for the 
linear power spectrum (3.13), it can be expressed as a weighted sum of the CDM+baryon 
and the neutrino density contrast auto- and cross-correlation power spectra, 

Pik, t) = f^Mk, r) + 2UUP,,,{k, t) + f^P^ik, r), (7.1) 

where 

P,bik, T)5D{k + k') = {5{k, T)5{k', r)), 
P,(fc, r)5D{k + k') = {S^ik, r)6''{k', r)), 

PcUk, r)6D{k + k') = {6ik, T)5^{k\ r)). (7.2) 

Assuming Gaussian initial conditions, these various power spectra, up to one-loop 
corrections, are given by [23] 

P^{k, r) = Pi{k, t) + [Pf )(A;, r) + 2P^'\k, r)], 

P,{k,T) = P!;{k,T), 

P^Uk, r) = Pt,M r) + P^l!\k, r), (7.3) 

taking the neutrino density perturbations to be linear at all times. 
The one-loop terms for the CDM-|-baryon auto-correlation are 

Pif (A;,r) = 2| d'q [Qf\k - q,q-r)f P^,{\k - qlT)P^,{q,r), (7.4a) 

P,f )(A;,r) = 3 P^,{k,T) J d'q Q^{'\k,q,-q;r)P^,{q,T), (7.46) 

where the prefactors "2" and "3" arise from summing all possible groupings of 
{6^{k)S^{k')) pairs to give {6^ {k)6^ {k')6^ {k")6^ {k'")) (e.g., [31]). The cross-correlation 
between CDM-|-baryons and neutrinos also contains a one-loop correction, 

P!:l'\k,T) = 3 P,i(A;,r) J d'q Qf\k,q, ~q;T)P^,{q,T), (7.5) 

M 3^ 

which is identical to the P^), {k,T) correction to the CDM+baryon auto-correlation 
except for the prefactor P^j^{k,T) instead of P^{k,T). This term is missing from 
Saito et a/.'s formulation, which assumed explicitly Pcbvik.r) = P^ij^{k,T) [21]. This 
assumption is not self-consistent, since a one- loop correction must be present in the 
CDM+baryon-neutrino cross-correlation, if the CDM+baryon density contrast has 
indeed been calculated to third order in perturbative expansion. 
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7.1. Explicit forms 
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Using equations (6.6) and (6.7) to evaluate Ql {k — q, q; r), and defining x = k-q/ (kq), 
r = q/k, and f] = y/l + r'^ — 2rx, we find for the "22" correction term (7.4a) 



2-Kk^ 
~W Jo 



drP^f^{kr, 



dxP^i,{kVl + r2 - 2rx, r) 



X 



[3Sir + 7(^2 + S3r^)x - lOS^rx 
(1 + — 2rxY 



212 



(7.6) 



with 



Si = ■^f{kv)Wii{kri,kr) - ay{kri,kr)f{kr)], 
5*2 = f{kr)[af^{kr], kr) + (yf.}{kri, kr)f{k'r])], 



af^{kr^M)[f{kr)-f{kii)l 
-af\kri,kr)f{kr), 



(7.7) 



and the functions /(g) and <J^^ {qi, ■ ■ ■ , qn) are defined in equations (3.2) and (6.3) 
respectively. In the f^—^O limit, we have 5'i_2,4 — ^ 1 and 5*3 ^ 0, which is the standard 
CDM-only result first given in references [32,33] and adopted in the analysis of [21]. 

Similarly, evaluating Q^i\k, q, —q; r) with the aid of equations (6.8) and (6.9), we 
obtain for the "13" terms (7.46) and (7.5) 

27r/t3 



63 



'Pcb,cbvi.k, t) 



drP^{kr,T) 







j dx 









(7.8) 



where 



and 



V = 
W - 

Hi 

H4 

h-- 

h-- 
h = 

h = 
L = 
M 



Qhr - l{l2 + hr'^)x + %hrx^ {lM-2R)r - {7Mr^-2R)x 



-a 



lOHirx - 7{H2 + H3r^)x^ + AH^rx^ 
^^{kM)[f{k) + f{kr)l 



7L 

f{kr)[af;!{k,kr) + af:^{k,kr)f{k)l 
- f{k)[af^{k,kr) + af^{kM)f{kr)l 

-laS{k,kr)f{k)f{kr), 



[7.9) 



-a. 



?^{k,kr)[f{k) + f{kr)], 



7 
■ —( 
6 

:f{kr)[a^^^{k,kr) + a!^^{k,kr)f{k)], 
■f{k)[a^^^{k,kr) + a^j{k,kr)f{kr)], 

\<yf^{k,kr)f{k)f{kr), 
ycrii^(/c, kr, kr)f{kr), 



7 



a'[f{k, kr, kr), 
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R = 9a^S{k,kr,kr)f{kr). (7.10) 

Again, these factors have been defined so that i^i, 2,3,4, -^1,2,3,4, and L, M, R tend to unity 
in the fu = hmit. Since these factors have no x-dependence, we can further simphfy 
equation (7.8) by performing the integration over x to obtain 

■12Ti 





'V 




1 ■ 


j dx 




vw 










~ 4 . 



I58T2 + lOOTgr^ - 42T4r' 



+ -(r2 - \)\T^r'' - n){mr^ + 2^) In 



1 + r 



(7.11) 



where 



T2 



2laf^{k, kr, kr)4l\k, kr)f{kr), 
42 
79" 



^f{kr) |al/(fc, kr, kr) [Aa\^{k, kr)f{kr) 

WS{k, kr)f{k)f{kr) - 34l\k, kr) 
+ 4f{k, kr, kr) [4? (A;, kr)[3f{k) + 9f{kr)] 

+Aa!Sik,kr)f{k)fikr)]}, 

T3 = § {-3ag(fc, kr, kr)a'S{k, kr)f{k)f{kr) 

+ af^{k, kr, kr) [4? (A;, kr)[?>f{kr) + 9/ (A;)] 

+ A[af^{k, kr) - af^k, kr)]f{k)f{kr) 
-Aa?j{k,kr)f{k)f{kr)]}, 

T, = 6a[f{k,kr,kr)aS^{k,kr)f{k), 
T, = f{k), 
n = f{kr), 

T7 = Qa^i{k, kr, kr)a!^i{k, kr), 

Ts = 21aSJ(fc, kr, kr)4S{k, kr)f{kr), (7.12) 

and Ti ... 8 — 1 as 0. The fu = hmit of equations (7.8) and (7.11) was first 

derived in [32,33] and then used in [21]. 

8. Results and discussions 

8.1. One-loop corrected power spectra 

Figure 1 shows various contributions to the total matter power spectrum at 2 = 1 
and 2; = 3 in a A mixed cold+hot dark matter (ACHDM) cosmology, assuming a 
total matter density fi^/i^ = 0.14, neutrino fraction f^, = 0.1, Hubble rate h = 0.72, 
and scalar spectral index Ug = 0.963. The total matter power spectrum has been 
normalised at the pivot scale ko = 0.002 Mpc~^ to the vanilla best-fit from the Wilkinson 
Microwave Anisotropy Probe five-year data, i.e., the best-fit amplitude of the curvature 
perturbations = 2.4x 10~^ [34]. Note that at the chosen redshifts the universe is very 
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k [h Mpc-i] k [h Mpc-i] 

Figure 1. Contributions to the total matter power spectrum (units: h'^ Mpc~^) 
at z = 1 (left) and z = 3 (right) for a ACHDM cosmology with ~ 0.1. 
Top: CDM+baryon auto-correlation Pc6(fc, t). Middle: CDM+baryon-ncutrino cross- 
correlation Pcbu{k,T). Bottom: Total matter power spectrum P(k,T). In all cases the 
linear contribution is shown in red/thin solid, the one- loop correction in green/thick 
solid (green/short dash when negative), and their sum in blue/long dash (blue/dotted 
when negative). All spectra have been divided by to facilitate comparison. The 
three vertical lines demarcate, from left to right, the upper limits in k for which the 
linear matter power spectrum is accurate to 1% or less, the one-loop corrected spectrum 
to 1% or less, and the one-loop corrected spectrum to 5% or less. 

nearly described by an Einstein-de Sitter model. However, the vacuum energy A term 
still has some residual effects, particularly on the logarithmic derivative of the linear 
growth function, f{k,T). It is therefore necessary to rescale f{k,T) so that f{k,r) 1 
as A; — > 0, in order for the expressions (7.6) to (7.8), (7.11) and (7.12) to be apphcable. 
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Observe that the P^^, {k,T) term derives from a positive definite integrand and is 
therefore always positive. The P'^^\k^T) term, on the other hand, is negative at most 
k values. The net effect of their sum on the CDM+baryon auto-correlation is that the 
power spectrum is first suppressed at /c ~ 0.01 — *■ 0.1 h Mpc~\ and then enhanced as 
we move beyond k ~ 0.1 h Mpc~^. 

Contrastingly, the cross-correlation spectrum between the CDM+baryon and 
neutrino components receives a one-loop correction only from P'^l\k,T). Like the 
P^^\k,T) term, P'^ll\k,T) is mostly negative. This causes the already free-streaming- 
suppressed power spectrum — recall the neutrino linear growth function falls off like 
at k ^ kps, see equation (3.12) — to be further suppressed, before turning negative at 
k r^0.2h Mpc"\ Beyond A; ~ 0.2 /i Mpc"\ however, the one-loop correction dominates 
over the linear term, resulting in an enhancement in the magnitude of Pcbu{k, t), but in 
the negative direction. 

The net correction to the total matter power spectrum follows essentially the 
same trend as the correction to Pcb{k,T), beginning with a small suppression at 
k < 0.1 h Mpc~^, and culminating in an enhancement at A; > 0.1 h Mpc~^. 

8.2. Regions of validity 

Perturbation theory is not expected to describe reality at very large k values, since 
any perturbative expansion must break down when the evolution of structures becomes 
fully nonlinear. A good rule of thumb is to take as the valid regime the range of k 
values at which the one-loop correction is smaller than the linear contribution [36]. 
More recent analyses show that for a ACDM cosmology, the one-loop corrected power 
spectrum deviates by less than 1% from A^-body simulation results provided that the 
dimensionless power spectrum, defined in this work as A^(A;,r) = 47r/c^P(A;, r), does 
not exceed ~ 0.4 [37]. Applying this criterion to the total matter power spectrum in 
our ACHDM scenario, we expect our one-loop corrected power spectrum to be accurate 
to better than 1% for A; < 0.2 h Mpc"^ at 2 = 1 and A; < 0.4 h Mpc"^ at ^ = 3. 
For completeness we also estimate a 5% accurate region following figure 2 of [37]: 
k <0.4:h Mpc"^ at z = 1, and k < 1 h Mpc^^ at ^ = 3. 

It is also interesting to ask up to what values of k is linear theory expected to 
provide an accurate description of the matter power spectrum. Comparing the total 
matter power spectrum computed from linear theory to that including the one-loop 
correction, we find 1% agreement between the two only at k < 0.09 h Mpc~^ for z = 1 
and k < 0.1 h Mpc~^ for z = 3. From this we conclude that the one-loop correction 
improves on linear theory to better than 1% accuracy in the k ranges 

0.09 < k/{h Mpc-i) < 0.2, z = 1, 

0.1 < A;/(/iMpc-i) <0.4, z = 3. (8.1) 
Figure 1 shows the various regions of validity discussed in this section. 




Figure 2. Relative differences between tlie total matter power spectra for a 
pure ACDM cosmology and three ACHDM models with = 0.1 (red/solid), 0.05 
(blue/dotted), and 0.01 (green/dash) at z = 1 (left) and z = 3 (right). Thick 
lines indicate results including the one-loop correction, while the linear results are 
represented by the thin lines. The three vertical lines indicate the maximum k values 
at which the linear and the one-loop corrected matter power spectra arc accurate to 
better than 1% and 5%. 



8. 3. Suppression due to neutrino free- streaming 

Figure 2 shows the suppression in the total matter power spectrum due to neutrino hot 
dark matter for three ACHDM models with f^, = 0.1, 0.05, 0.01, relative to the case with 
fu = 0, i.e., 

P{k,a) Pf^=o{k,a) 

Observe that the relative decrease of small scale power due to neutrino free-streaming 
in general exceeds the amount of suppression predicted by linear theory, once nonlinear 
corrections have been included. Within the 5% accurate region, it is clear that the 
suppression easily exceeds the canonical linear suppression factor of ~ 8/,^.§ This 
enhanced suppression has been observed in multi-component A^-body simulations [18], 
which support an asymptotic suppression factor of ~ 9.8f^. Furthermore, compared to 
the linear results, the one-loop correction leads to a small increase in relative power at 
just below k ^ 0.1 h Mpc~^, before the enhanced suppression sets in, a feature that 
seems to be present also in figure 4 of [18] 

A more rigorous comparison between our perturbation theory results and the A^- 
body results of [18] is not possible at present, since our results are valid only for 2 > 1, 
while reference [18] gives their results at 2; = 0. 

§ The linear suppression factor reaches a maximum of ~ 8/,^ only for small values of f^. For, e.g., 
fiy = 0.2, the linear suppression factor is ^ 4.5/,^ [35]. 
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Figure 3. One-loop corrections to the CDM+baryon auto-correlation (top) and the 
CDM-f baryon-neutrino cross-correlation (middle) at z = 1 (left) and z = 3 (right) for 
a ACHDM cosmology with ~ 0.1. In the top and middle rows, the red/thin solid 
lines indicate the linear contribution, while the green/thick solid lines denote the one- 
loop correction (green/short dash when negative). The blue/long dash lines represent 
an approximation to the one- loop correction (blue/dotted when negative), computed 
by setting 51.2.4 ~ Ti ...^^ = 1 said S3 = 0. The fractional error incurred in the total 
matter power spectrum by this approximation is shown in the bottom row in red/solid. 
Also plotted in the bottom row in blue/long dash is the fractional error incurred by 
neglecting the one-loop correction to the CDM+baryon-neutrino cross-correlation. 



8.4- Further approximations? 

There are two sources of deviation from the standard CDM-only case in the one-loop 
correction terms (7.6), (7.8) and (7.11). The first is encapsulated in the linear power 
spectra P^{k,T) and P^j^{k,T)] the second in the factors 5'i_..._4 and Ti^...^8 defined in 




Figure 4. Same as figure 3, but for = 0.05. 



equations (7.7) and (7.12), which depend on the hnear growth functions. The hnear 
power spectra are readily calculable with a sophisticated Boltzmann code such as CAMB. 
With some extra but minor effort the same is true also for the factors 5*1^. ..^4 and Ti ... §. 
Nonetheless, it is tempting to simply assume 5'i,2,4 = Ti^...^s = 1 and 5*3 = 0, as was done 
in [21]. In this section we examine the validity of this approximation. 

Figures 3, 4 and 5 show the correction terms Pcb^"* + ^chu computed 

under the assumption of 51,2,4 = Ti,...^^ = 1 and 53 = for three ACHDM models with 
= 0.1,0.05,0.01 respectively. Compared with their correct forms, we find deviations 
as large as a factor of twenty, sometimes even accompanied by a sign flip. However, 
such large deviations are confined to a region at k 0.1 h Mpc~^, where the one- 
loop correction terms are in any case subdominant to the linear contribution; the net 
fractional error incurred in the total matter power spectrum turns out to never exceed 




Figure 5. Same as figure 3, but for = 0.01. 



the 1% level for our choices of /j, and redshifts. Future cosmological probes will generally 
require an accuracy of ~ 1% in the matter power spectrum in order not to bias parameter 
estimation. On this basis, we conclude that setting 5'i^2,4 = ^i,- -,8 = 1 and 6*3 = is a 
tolerable simplification. 

Also shown in figures 3, 4 and 5 are the consequences of dropping the P^l^\k,T) 
correction to the CDM+baryon-neutrino cross-correlation. As expected, the importance 
of this term increases as we increase For f,y = 0.1, the maximum contribution of 
P^l\k,T) to the z = 1 total matter power spectrum is 2.5% at A; ~ 0.4 h Mpc~"^. For 
/(, = 0.05, the contribution is essentially halved. Thus for fi, < 0.05, dropping the 
P'^l\k, t) term will unlikely bias our results. Nonetheless, P'^l\k, r) can be computed 
at essentially no extra expense to the user because of its similarity to the (mandatory) 
P^jj {k,T) contribution, cf. equation (7.8). Hence there is no real need to resort to 
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approximations in this instance. 

9. Conclusions 

In this paper we have presented the first rigorous and systematic derivation of the 
one-loop correction to the large scale matter power spectrum in a mixed dark matter 
cosmology with subdominant massive neutrino hot dark matter. 

Beginning with the relevant equations of motion, we find that by invoking an 
"adiabatic" approximation, accurate to better than ~ O.l/j,, higher order corrections to 
the CDM+baryon density contrast and velocity field can be rendered into a form nearly 
identical to that for a pure-CDM cosmology. The interaction kernels and their recursion 
relations also exhibit striking similarities to their standard CDM-only counterparts, but 
contain additional dependences on the neutrino energy density fraction fi, and the linear 
growth functions of the incoming wavevectors. These results, generalised to nth order 
in perturbative expansion, are summarised in equations (6.1) to (6.4). 

Using these approximate solutions we compute the usual "22" and "13" one-loop 
correction terms to the matter power spectrum. As in the standard CDM-only case, 
these correction terms take the form of integrals over the wavevector q of the linear 
power spectrum P^{q,T) multiplied by the interaction kernels. In addition to the 
corrections to the CDM-|-baryon auto-correlation, we also find a one-loop correction 
term for the cross-correlation between the CDM+baryon and the neutrino components 
which was previously neglected. These correction terms appear in their evaluated and 
most simplified form in equations (7.6) to (7.8), (7.11) and (7.12). 

Evaluating these expressions numerically, we find that nonlinear corrections to the 
large scale matter power spectrum can enhance the suppression of small scale power 
due to neutrino free-streaming relative to the = case to beyond the canonical 
linear suppression factor of ~ 8/^. This enhanced suppression has been observed in 
multi-component A^-body simulations [18]. 

As said, the interaction kernels contain hitherto unaccounted dependences on 
and the linear growth functions. Neglecting these dependences in principle generates 
large deviations in the one-loop corrections. However, since these deviations occur at 
wavenumbers at which the linear contribution dominates over the correction terms, their 
net effect on the total matter power spectrum never exceeds 1%. Future cosmological 
probes will require an accuracy of ~ 1% in the matter power spectrum in order not to 
bias parameter estimation. We have thus verified the validity of the approach of [21]. 

An important assumption in our present treatment is that the neutrino density 
perturbations have been taken to remain linear at all times. Although for realistic 
values of fi, this can be justified by A^-body simulation results [18], a truly complete 
analysis of higher order corrections to the clustering statistics of the large scale structure 
distribution in the presence of massive neutrinos should include also a proper account 
of nonlinear neutrino evolution. We defer this investigation to a future publication. 

Finally, as noted in reference [21], although higher order perturbation theory 
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appears at first glance to have a limited range of validity — we expect our one-loop 
corrections to improve on linear theory to better than 1% accuracy in the region 
0.1 ^ k/(h Mpc~^) < 0.4 at 2; = 3, it does enable an approximate factor of four increase 
in the maximum usable wavenumber in a data set. This is equivalent to a factor of 64 
gain in the number of independent Fourier modes, or an eight-fold gain in statistical 
power for a fixed survey volume. Such an improvement is no small feat, and may very 
well be just what we need to detect neutrino dark matter. 
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